/*
 * EBFV1-hyperbolic.cpp
 *
 *  Created on: 09/01/2009,2012
 *      Author: rogerio
 */

#include "EBFV1_hyperbolic_adaptative.h"
#include "exportVTK.h"

namespace PRS{
	EBFV1_hyperbolic_adaptative::EBFV1_hyperbolic_adaptative(){}
	EBFV1_hyperbolic_adaptative::EBFV1_hyperbolic_adaptative(PhysicPropData* _pPPData, SimulatorParameters* _pSimPar, MeshData* _pMData, OilProductionManagement* _pOP, TMesh* _pTMesh_data){
		pSimPar = _pSimPar;
		pOPManager = _pOP;
		pPPData = _pPPData;
		pHOApproximation = new HighOrderApproximation(pSimPar);
		initial_vel_ts_old = true;
		print_ATE_Data = true;
		vel_ts_counter = 0;
		sat_ts_counter = 0;
		DVTOL = pSimPar->getDVTOL();
		createDimensionLessFactors();
	}

	EBFV1_hyperbolic_adaptative::~EBFV1_hyperbolic_adaptative(){
	}

	double EBFV1_hyperbolic_adaptative::solver(TMesh* pTMesh_data, double &timeStep){
		PetscPrintf(PETSC_COMM_WORLD,"Adaptative Hyperbolic solver...\n");
		double startt = MPI_Wtime();
		int dim = pTMesh_data->getMeshDim();
		int ndom = pTMesh_data->getNumDomains();
		double sat_ts_sum = .0;

		resetVariables();
		while ( allowAdaptativeTimeAdvance ){
			// reset logical variables which control programming flux.
			this->resetTimeStep();

			// For High-order approximations
			if ( pSimPar->useHOApproximation() ){
				calculateSaturationGradient(pTMesh_data);
				pHOApproximation->getNodeSL_Ptr()->defineSlopeLimiters();
			}

			//	calculateSaturationGradient();
			for (int dom=0; dom<ndom; dom++){
				if (allowVelocityCalculation){
					calculateVelocityField(dom,dim,pTMesh_data);
				}
				resetNodalNonviscTerms(dom,pTMesh_data);
				calculateIntegralAdvectiveTerm(dom,dim,pTMesh_data);
			}
			//pMData->unifyScalarsOnMeshNodes(PhysicPropData::getNonViscTerm,PhysicPropData::setNonViscTerm,pGCData,0);

			sat_ts = getTimeStep();														// take minimum time-step calculated by domain
			pSimPar->correctTimeStep(sat_ts);											// correct time-step value to print out the desired simulation moment
			pSimPar->setAccumulatedSimulationTime(sat_ts);								// AccSimTime = AccSimTime + timeStep
			vel_ts = calculateNewImplicitTS(pTMesh_data);								// calculate implicit time-step
			saturationTimeAdvance(pTMesh_data,sat_ts);									// Calculate saturation field: Sw(n+1)
			pSimPar->printOutVTK(pTMesh_data,exportSolutionToVTK);						// Output data (VTK)
			if (pSimPar->rankHasProductionWell()){										// The following lines below will be condensed to a function member call and it will belong to EBFV1_hyperbolic
				pOPManager->printOilProduction(sat_ts,pSimPar->getAccumulatedSimulationTime(),pSimPar->getSimTime());
			}
			sat_ts_sum += sat_ts;
			sat_ts_counter++;
			if (sat_ts<.0){
				throw Exception(__LINE__,__FILE__,"Negative time-step\n");
			}
		}

		vel_ts_counter++;
		timeStep = sat_ts_sum;

		if (!P_pid()){
			printAdatptativeTimeEvaluationData();
			cout << "done.\n";
		}
		double endt = MPI_Wtime();
		return endt-startt;
	}

	// Compute a new implicit time-step (vel_ts) if necessary. This happens either at the beginning of simulation or when sum of sat_ts is greater than DT.
	double EBFV1_hyperbolic_adaptative::calculateNewImplicitTS(TMesh* pTMesh_data){
		PrintForDebugging("calculateNewImplicitTS ...\t");

		if (allowVelocityCalculation){
			// stop velocity field be evaluated during saturation advance.
			allowVelocityCalculation = false;

			// set vel_ts_old as sat_ts just once at the beginning of simulation.
			if (initial_vel_ts_old){
				vel_ts_old         = sat_ts;
				initial_vel_ts_old = false;
			}

			// make implicit vel_ts_old variable dimensionless
			vel_ts_old *= time_factor;

			// implicitTS returns a dimensionless delta T
			dv_norm = calculateVelocityVariationNorm(pTMesh_data);

			// Update new implicit time step: vel_ts
			if (fabs(dv_norm) < 1e-10){
				cout << "Warning: division by zero (dv_norm = " << dv_norm << ")\n";
			}
			vel_ts = (DVTOL/dv_norm)*vel_ts_old;
			//PetscPrintf(PETSC_COMM_WORLD,"DVTOL: %f dv_norm: %f  DVTOL/dv_norm: %f time_factor: %f\tvel_ts_old: %f\tvel_ts: %f\n",DVTOL,dv_norm,DVTOL/dv_norm,time_factor,vel_ts_old,vel_ts);

			// dimensionless DT ratio
			double RT = (double)(vel_ts/vel_ts_old);
			if (RT>1.25){
				vel_ts = 1.25*vel_ts_old;
			}
			else if ( RT<0.75 ){
				vel_ts = 0.75*vel_ts_old;
			}

			// gives implicit vel_ts variable dimension of time
			vel_ts /= time_factor;
			vel_ts_old = vel_ts;

			// It doesn't make sense vel_ts be less than sat_ts
			if (vel_ts < sat_ts){
				vel_ts = sat_ts;
				allowAdaptativeTimeAdvance = false;
			}
		}
		// Verify explicit advance compared to implicit advance: sum(sat_ts)<=vel_ts
		cumulativeExplicitTS += sat_ts;				// summation of all sat_ts within the implicit time-step (vel_ts)
		verifyExplicitAdvanceForImplicitTS();

     	// Verify explicit advance compared to simulation time. AccSimTime = AccSimTime + timeStep. Summation of all sat_ts for the whole simulation
		///pSimPar->setAccumulatedSimulationTime(sat_ts);
		verifyExplicitAdvanceForSimulationTime();
		PetscPrintf(PETSC_COMM_WORLD," DONE\n");
		return vel_ts;
	}

	// Compute a new implicit time-step as function of the last implicit time-step, velocity norm and DVTOL. DVTOL comes from the input data file 'numeric.dat'.
	double EBFV1_hyperbolic_adaptative::calculateVelocityVariationNorm(TMesh* pTMesh_data){
		PrintForDebugging("calculateVelocityVariationNorm ...\t");

		const int dim = pTMesh_data->getMeshDim();
		const int ndom = pTMesh_data->getNumDomains();
		const int numGEdges = pTMesh_data->getNumGEdges();
		double v_var[3];
		double Dv_dom[ndom];

		int k,i,j,dom,nedges,n;
		// loop over domains
		k = 0;
		for (dom=0; dom<ndom; dom++){
			Dv_dom[k] = .0;
			pTMesh_data->edge_list[dom].getSize(nedges,n);
			for (i=0; i<nedges; i++){
				EdgeData* edata = pTMesh_data->edge_list[dom].getValue(i,0);				// get next domain's edge
				// make dimensionless velocity
				for (j=0; j<dim; j++){
					v_var[j] = (edata->vel_new[j] - edata->vel_old[j]);
				}
				double innerp = dot(v_var,edata->Cij,dim);
				Dv_dom[k] += pow(fabs(innerp)/edata->Cij_norm,2);
			}
			//Dv_dom[k] = P_getSumDbl(Dv_dom[k]);
			Dv_dom[k] /= (double)numGEdges;
			k++;
		}
		// take an average of velocity norm for all domains
		double aux = (double)sum(Dv_dom,ndom)/ndom;
		double vnorm = sqrt( aux );
		PrintForDebugging("Done.\n");
		return vnorm;
	}

	void EBFV1_hyperbolic_adaptative::resetVariables(){
		cumulativeExplicitTS = .0;
		allowVelocityCalculation = true;
		allowAdaptativeTimeAdvance = true;
	}

	void EBFV1_hyperbolic_adaptative::verifyExplicitAdvanceForImplicitTS(){
		if ( cumulativeExplicitTS >= vel_ts ){
			sat_ts = sat_ts - (cumulativeExplicitTS - vel_ts);
			allowAdaptativeTimeAdvance = false;
		}
	}

	void EBFV1_hyperbolic_adaptative::verifyExplicitAdvanceForSimulationTime(){
		double accSimTime = pSimPar->getAccumulatedSimulationTime();
		double simTime = pSimPar->getSimTime();
		if ( accSimTime >= simTime ){
			sat_ts = sat_ts - (accSimTime - simTime);	// correct sat_ts (last step)
			allowAdaptativeTimeAdvance = false;
			pSimPar->stopSimulation();
		}
	}

	void EBFV1_hyperbolic_adaptative::printAdatptativeTimeEvaluationData(){
		if (print_ATE_Data){
			int size = P_size();									// number fo processes
			double dvtol = pSimPar->getDVTOL();			// dvtol value
			string path = pSimPar->getOutputPathName();	// where file will be placed
			char fname[256]; sprintf(fname,"%s_AdaptTimeData-DVTOL__%.2E__%dp.csv",path.c_str(),dvtol,size);
			fid.open(fname);
			if ( !fid.is_open() ){
				char msg[256]; sprintf(msg,"%s could not be opened or it does not exist.\n",fname);
				throw Exception(__LINE__,__FILE__,msg);
			}
			fid << "Time-step ImpExp_TimeStep-ratio velocity_norm\n";
			print_ATE_Data = false;
		}

		if (!P_pid()) {
			double accSimTime = pSimPar->getAccumulatedSimulationTime();
			double simTime = pSimPar->getSimTime();
			double time_ratio = (double)sat_ts_counter/vel_ts_counter;
			// convert numeric data into a string
			char lineStr[256];
			sprintf(lineStr,"%f %f %f",accSimTime/simTime,time_ratio,dv_norm);
			string theString(lineStr);
			replaceAllOccurencesOnString(theString,1,".",",");
			fid <<  theString << std::endl;
		}
	}

	void EBFV1_hyperbolic_adaptative::createDimensionLessFactors(){
		double L = .0;
		double H = .0;
		pSimPar->getReservoirGeometricDimensions(L,H);
		double Q = pSimPar->getTotalInjectionFlowRate();
		time_factor = Q/(pSimPar->getPorosity(3300)*pow(L,2)*H);
		vel_factor = L*H/Q;
	}
}
